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1 Introduction 

Modern perturbative computations typically involve large numbers of Feynman diagrams 
and/or Feynman integrals, whose successful evaluation is greatly facilitated by the use of 
computer algebra. Indeed, in the field of multi-loop calculations (in zero-temperature field 
theories) the algorithmic formulation is in a quite mature state, as witnessed by numerous 
higher-order results, such as e.g. 5-loop results for the QED beta function [1], complete 3- 
loop and parts of 4- loop terms in the electron anomalous magnetic moment [2], or 4- loop 
contributions to electroweak precision observables [3], just to name a few. 

Such calculations are typically performed in a sequence of steps: (a) generation of the 
complete set of diagrams and counter-terms contributing to the observable under study; (b) 
application of Feynman rules, necessary projectors and traces, performing Lorentz algebra and 
scalarization; (c) mapping onto a set of integral templates; (d) reduction to a few ("master") 
integrals; (e) expansion in epsilon. 

For each of these distinct steps, a variety of systematic algorithmic methods and tools 
can be applied, the most commonly used ones being: (a) graph theory (efficiently coded in 
the package QGRAF [4]); (b,c) computer-algebra systems (such as e.g. FORM [5], Ginac [6], 
Mathematica [7]); (d) integration- by-parts methods (pioneered by Chetyrkin and Tkachov 
[8], formalized by Laporta in [9] and implemented e.g. in the public packages Air [10], FIRE 
[11], Reduze [12]); (e) diff'erence equations [9], diff'erential equations (see, e.g. [13]), harmonic 
polylogarithms and -sums [14] (where the expansion can be automatized for simple 0-scale 
problems [15], but is non-trivial for multiple-scale problems), graph polynomials [16], sector 
decomposition [17]. 
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For problems in finite temperature field theories, however, a similar computer-algebraic 
approach to higher-order calculations is much less developed, although a limited number of 
results up to the three-loop [18-20] and even the four-loop level [21] do exist. Whenever 
gauge theories are treated, owing to the structure of gauge-field propagators and -vertices, 
most authors choose to work in a fixed (typically Feynman) gauge, in order to reduce the 
complexity of the calculation. 

While most of these works were performed in a more traditional way, there is no major 
obstacle in carrying over the systematics of the modern computer-algebraic developments 
mentioned above to finite temperature calculations, which as a consequence abolishes the 
need to work in fixed gauges. In fact, for the steps (a) to (d) - i.e. diagram generation / 
algebraic simplifications / mapping / reduction - the methods from zero-temperature field 
theory can be generalized directly, with a few minor modifications. It is only step (e) - i.e. 
the epsilon expansion of master integrals - that resists full automation, owing to the fact 
that the available techniques for sum-integrals which are needed in the thermal setting are 
much more limited than those for pure continuum integrals (for which at least numerical 
methods are guaranteed to work). For first attempts in treating whole classes (as opposed to 
treating them one-by-one as has been the state of the art previously [18, 21, 22]) of non-trivial 
sum- integrals, see [23]. 

In this paper, we want to demonstrate the utility of computer-algebra methods for thermal 
field theories, which mainly concerns the integral reduction step (d). Concentrating on the 
concrete example of the 3-loop free energy of hot QCD, we briefly introduce the corresponding 
observable in Sec. 2 and then explain the specifics of the (IBP) reduction method generalized 
to finite temperature in Sec. 3. Sec. 4 contains the result of our diagrammatic calculation, 
confirming the known 3-loop result [19], however in general covariant gauges. In Sec. 5 we 
conclude. The Appendix lists some useful sum-integrals that are needed for the final result. 



2 Thermodynamic observables 

QCD equilibrium properties, such as its free energy density F, are encoded in the logarithm 
of the partition function 

F = In jv[A^,i^,i^] exp fj^ Jd'^xCQCB^ . (2.1) 

Here, V = fd'^x denotes the spatial volume occupied by the system {d = 3 — 2e), and we 
work in the imaginary time formalism, defined on a ((i-l-l)-dimensional Euclidean space with 
a compact temporal coordinate r with period 1/T set by the temperature. Bosonic/fermionic 
(gluon/quark) fields are periodic/antiperiodic functions of r. The (Euclidean) QCD La- 
grangian reads 

1 _ 

^QCD = ^ F^,F^, + 5Z + "^^ + ^'"*) ^ ' (2.2) 

i=l 
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with field strength tensor F^^^ = dfj_A'^ — dyA'^ + g f"''"^ A^^Af, and covariant derivative Z?^ = 
+ igA'^T"' . We will set quark masses rrii and chemical potentials to zero here, and work 

with the gauge group SU(A^c), where Ca = A^c, TY(r"r'') = \ 6"^ and T^T" = Cpl = 1- 
It turns out that the weak-coupling expansion of Eq. (2.1) is nonanalytic in the strong 
coupling constant Os = g"^ /Air. The physical reason is that, due to multiple interactions in 
the thermal medium, screening masses are dynamically generated for all massless particles, 
resulting in a multi-scale system. In fact, at high temperatures, asymptotic freedom guaran- 
tees a small gauge coupling g. In this regime, QCD develops a hierarchy of three momentum 
scales vrT ^ gT ^ g'^T /t:, whose effect can be most transparently accounted for in an effec- 
tive theory setup [24]. Systematically integrating out the largest ("hard") scale vrT and the 
second-largest ( "soft" ) scale gT in turn, one obtains a dimensionally reduced effective theory 
[25] for the smallest ( "ultrasoft" ) scale, which has been dubbed magnetostatic QCD (MQCD) 
[19]. The leading term of the MQCD action turns out to be a 3-dimensional pure Yang-Mills 
theory, which is confining and therefore has to be treated with suitable non-perturbative 
methods [26]. 

Utilizing the effective theory setup, the QCD pressure can be expressed as [19] 

PQCD(r) = - ,lim F = Phard(7') + Psoit{T) + Pultrasoft (?") , (2-3) 

where phard and psoft are perturbatively computable matching coefficients which account for 
contributions from hard and soft momentum scales, respectively. Defining these perturbative 
matching coefficients in the MS scheme, the remaining contribution 

J3uitra.oft(r) = j^lim^ ^ In jv[A,] exp jd'^x Cmqcb^ | (2.4) 

entails the effective Lagrangian £mqcd = jF^j^ij + where F^j = diA'j — djAf + 
guf'^^'^A^Aj contains the 3-dimensional gauge coupling gu- In contrast to the dimensionless 4- 
dimensional gauge coupling g, the dimensionality of g^^ is GeV, such that for dimensional rea- 
sons Puitrasoft ~ 2^ 5m- ^he other hand, perturbative matching yields 5^ = g^T (1 + 0{g)), 
so Puitrasoft plays a role in pqcD starting at 0{g^) only, which is beyond the precision needed 
for our investigation (see, however [26]). 

The two matching coefficients in Eq. (2.3) are defined in analogy to Eq. (2.4), but with ac- 
tions jQ^^dr Jd'^xCqcD and jd'^x Ceqcd for phard andpsoftj respectively. The latter depends 
on the so-called 3d electrostatic QCD (EQCD) Lagrangian containing a massless gauge field 
Ai and a massive adjoint scalar Aq, and whose structure results from integrating out the hard 
scales from QCD, yielding the leading terms >Ceqcd = j ^ij^ij + Tr[A^o] [A^o] + m|TrAg + 

(TrAg)^ -|- A^^Tr^Q + . . . . For a review on the status of the different contributions we 
refer to [27]. 

It turns out that at present, the bottleneck is the evaluation of higher orders in Phard(^)) 
which can be obtained by adding all vacuum diagrams in 4-dimensional thermal QCD, eval- 
uated in the "naive" perturbative sense, i.e. by regularizing ultraviolet as well as infrared 
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divergences in dimensionally (it is the effective theory setup Eq. (2.3) that properly accounts 
for infrared effects [28] that need to be resummed). The need to ultimately perform 4-loop 
computations within this setting is our main motivation to proceed with computer-algebraic 
methods as far as possible. As mentioned in the introduction, this concerns mainly the task 
of integral reduction, which we will now discuss, and then apply to the problem of 3-loop 
corrections to Phard(^)) enabling a comparison with known results from the literature. Going 
a step beyond the Feynman-gauge treatment of [19], we work in covariant gauges with gluon 
propagator 

(2.5) 

explicitly keeping the gauge parameter ^ (note that in our convention, ^ = 0(1) corresponds 
to Feynman (Landau) gauge, respectively), and demonstrating its cancellation in the sum of 
diagrams. 



ab 



?i P P 

p2 ~ ^Jp2y2 



3 Integral reduction 

For the sake of concreteness, and to avoid too generic notation and proliferation of indices, 
let us take the problem of three-loop sum-integral reduction as an example here. While it 
should be understood that the methods introduced below are independent of this choice, let 
us note that this is in fact the first non-trivial loop order in the case of sum-integrals (cf. 
Appendix A), and also precisely the level of the computation displayed in Sec. 4. 
In general, 3-loop vacuum-type sum-integrals can be written in the form 

.a/37 =4 (^o)" {Q^f {Ror .... 

abcdef;c,c,c. - ^^^^ [p2]a [Q2]b [p2]c [(p _ g)2]d [(p _ p)2]e [(g _ p)2]/ ' ^ ' ) 

where P^ = {Pq)^ + = ([2np -|- Cp]7rT)^ -|- are bosonic (fermionic) loop momenta for 
Cj = (1), and where the indices a . . . f € Z and a ... 7 € Nq. The sum-integral symbol in 
Eq. (3.1) is a shorthand for 

d'^p 



Po 

where // is the minimal subtraction (MS) scheme scale parameter, we take d = 3 — 2e, and the 
sum is over all integers Up £ 'Z. Hence, the set of indices of the sum- integral / enumerates all 
possible structures that can occur in a particular (3-loop) computation of finite-temperature 
Feynman integrals without further external momentum scales (a generalization to n-point 
functions or to internal masses and/or chemical potentials is straightforward, but let us focus 
on the problem at hand here). 

Generic integration- by-parts (IBP) relations then provide linear relations between the set 
of sum-integrals / of Eq. (3.1), using that the integral of a total derivative (here with respect 
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to the spatial loop momenta only) vanishes in dimensional regularization, 



= 



t dpJ,ip„Po,...)g{P,...) 



(3.3) 



Jp... 



with arbitrary function fi (and where g denotes an integrand of the type Eq. (3.1)). The 
linear relations are obtained by choosing the function fi, working out the derivatives, and re- 
expressing the result in terms of the generic form Eq. (3.1). From this point on, it is clear that 
the well-established (zero-temperature) algorithms that systematically solve such systems of 
linear IBP relations can be taken over. In practice, providing a unique ordering relation among 
the sum-integrals /, we use two Laporta-type [9] algorithms (one programmed in FORM [5], 
as well as one in Ruby [29]; the latter code is used as a cross-check on the former, setting the 
dimension d to a numerical value, utilizing the speed of numerical Gaussian elimination, and 
avoiding polynomial algebra). 

Furthermore, to reduce the number of relations, it is useful to exploit symmetries among 
the I, which can be generated by linear shifts of loop-momenta. For example, due to the sym- 
metry of the generic sum-integral in Eq. (3.1) it is sufficient to consider, from all 2'^ = 8 possi- 
bilities, the three cases {cp, Cq, c,.) = {(000), (100), (110)}. This mapping can be automated by 
systematically treating linear relations originating from momentum shifts on the same footing 
as the IBP relations. A typical relation originating from such shifts is /oooiii ooi ~ -^iiiooo iio- 
In general, however, for non-zero coefficients a, /?, 7 and/or negative values among the indices 
a — / (i.e. scalar products in the numerator), these relations have more than one term on the 
right-hand side. 

One could suspect that there exist additional IBP-type relations that take into account 
the structure of the Matsubara sums, which Eq. (3.3) has not sampled yet. This is not the 
case, however, as we will show here. In fact, acting with the operator TOt on a generic L-loop 
sum-integral (assuming massless propagators, which however can be of bosonic or fermionic 
type), one can either use the fact that, due to the absence of any other dimensional scale, 
the dimension of the sum-integral is carried by the scale T only, or apply the derivative to 
the explicit (in front of each sum) and implicit (in the Pq = {2n + cJttT of the integrand) 
occurrences of T directly 



It turns out that Eq. (3.4) carries the same information as the sum of the "diagonal" IBP 
relations 




Pio Plo 




(3.4) 




g{Pi,...,PL) 



(3.5) 
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and therefore does not need to be considered. We have checked this exphcitly up to four 
loops, for vacuum sum- integrals. 

One can however use an additional set of relations that derive from the sum-part of 
the sum-integrals. These additional relations essentially mix bosonic with fermionic sum- 
integrals, and are based on scaling arguments, such as used e.g. in Appendix B of [30] (see 
also [31]). One first re-scales the spatial integration momenta of a given bosonic integral as 
Pi — >■ 2pi and then partitions the Matsubara sums as 

E = E+E- (3-6) 

nGS even odd 

In practice, this provides a few linear relations among different bosonic and fermionic master 
sum- integrals, which remain after systematic use of the IBP relations Eq. (3.3). The simplest 
example is the relation between 1-loop tadpoles, as shown in Eq. (A. 2). For the 3-loop 
calculation presented in Sec. 4 below, we obtain (letting /nooii abc = ^abc for brevity) 

1 

^000 



PQR P'Q' {P-RY{Q-RY 
V V V ^ 1 



^ ^ ^ ioar [(2n„^r)2 + p2] [(2n„^r)2 + q^] . . . 

2^ J [(n„7rT)2 + p2] [(n„7rr)2 + q2] 



l3d 



-^000 + -^001 + -^010 + hm + hw + -^loi + -^oii + -^iii 

= 23'^-^ (/ooo + 6 /ool + /no) , (3.7) 

where in the third line we have scaled the spatial momenta, in the fourth line considered all 
cases of Eq. (3.6) (cubed), and finally exploited symmetries of this simple basketball-type 
sum-integral. Altogether, we therefore have the linear relation (anticipating the notation of 
Eqs. (4.11)ff for master integrals) 

= (l-28-3'^)/ooo + 6/001 + /no = {l-2''-^'^)B2 + &F^ + F-r (3.8) 

between one bosonic and two fermionic integrals, which will allow us to reduce the basis of 
master integrals by one. We feed these types of linear relations into our IBP system as well. 

4 Results of diagrammatic calculation 

Let us now evaluate the coefficient Phard(/')i as defined in Sec. 1, to three- loop order. In 
doing so, we work in covariant gauges (cf Eq. (2.5)), aiming at proving gauge parameter 
independence as well as confirming the corresponding Feynman-gauge result of [19]. Working 
in dimensional regularization with d = 3 — 2e, let us rewrite bare quantities as 
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Figure 1. "One-loop" contributions to the QCD pressure. Wiggly/dottcd/solid lines denote glu- 
ons/ghosts/quarks, respectively. The alphabetical index on the prefactor labels the diagram. 




Figure 2. Two-loop contributions to the pressure. Notation as in Fig. (1) 

The leading ("one-loop" - or, more precisely, the logarithms of Gaussian path integrals 
over the quadratic parts of the action, given by the sum over logarithms of momentum-space 
eigenvalues of the corresponding matrix kernels) terms of Fig. (1) can be expressed in terms 
of basic logarithmic integrals, which in turn can be related to the more conventional one-loop 
sum-integrals over propagator structures given in Eq. (A.l). For the bosonic case. 

Here, we noticed that / depends on the scale T only, such that on the one hand its derivative 
gives the overall dimension, while on the other hand the derivative can be applied directly 
(hitting the explicit T in the sum-integral as well as the implicit one in P^). In complete 
analogy, one gets for the fermionic case 

i ^ HP') = \ n = I (2-'^ - 1) ii , (4.3) 

such that the individual diagrams of Fig. (1) contribute as 

p':' = -'-^d,ll '='T^%d,, (4.4) 

p'^ = \d^ll '=-^^^rfA, (4.5) 

^1 = \ (2-^-1) N,C^ll '=^' T^I^ N,Ca , (4.6) 

whose sum gives the well-known (QCD version of the) Stefan-Boltzmann law. One can 
see clearly the effect of the ghosts here, which cancel half of the result of the pure gluonic 
contribution. 

After reduction of the 2-loop diagrams of Fig. (2), we obtain d-dimensional bare results, 
expressed in terms of bosonic as well as fermionic 1-loop tadpoles (as has been pointed out 
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Figure 3. Three- loop diagrams contributing to C^- 



e.g. in [19], all 2-loop sum-integrals factor into products of two 1-loop cases; this obser- 
vation we reproduce via IBP, for all possible 2-loop vacuum sum-integrals, also of different 
dimensionality needed for other computations; see also Appendix A) 

d{d + iy 



[d] 



Ph 



g^dxCx 
g^dANf- 



d{d-3) ^2 



2 ^ 



tOtO 

-'l-'l ) 



1 



d 



16 
1 

' 4 
2/V 



rO tO 



g'dANr 



d-1 



3 ) /?/? 



(4.7) 
(4.8) 
(4.9) 
(4.10) 



Turning now to the three-loop case, the IBP reduction leaves us with a set of bosonic 
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-D3 = 



000 
12000;000 - 

000 

21111-2;000 
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-'l -'l -'2 ) 

B4 = i 
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110012:000 5 



B2 = T 
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000 

110011;000 ' 



7-000 

'11-2013:000 ' 



(4.11) 
(4.12) 



as well as fermionic master sum- integrals of the form Eq. (3.1) (note however the linear 
relation Eq. (3.8) that can still be used to eliminate one of B2,Fq,Fi) 



Fi 


— 7-000 

= -'l21000;001 


= /? 


1^2 




(22- 


-d _ 


l)Bi , 


F2 


_ 7-000 

= -'ll2000;001 


= /? 




11 = 


(2^- 


-d _ 
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F3 


_ 7-000 

= -'ll2000;011 
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-d _ 
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_ 7-000 
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(4.13) 
(4.14) 
(4.15) 
(4.16) 
(4.17) 



In terms of these masters, we obtain d-dimensional bare results from the three-loop diagrams 
depicted in Fig. (3): 



\h+i+j+k] 
Ph. 



d{2 + 33(f - 12d2 + d^) ^4 ^ d{-6 - 13d + 3^^) 



384 



32 



- 8 - 



2 + 41d + 56d^ - 8d^ ^ 13 + 115d + 38^2 - 50^^ + 3 - lU + A(f \ 

^ 64 ^ 8(d-3)((i-5) ^ 8 

/ 22 - 13d + 3d^ ^2 11 - 19d + ^ 5 - lOd + A(P \ 

^ ^\ 768 ^ 192(d - 3) 32 

/ (d-l)(3(i-10) ^2 , 2S-nd + M\ , 1 

-^^V ^ + 96(d-3) ^ + I6j + 

d - 3 ,2 1 A D 1 ^2 1 



64 ' V 768 192 ' 32 



W 4^ ^2/^o ^ '^(rf-S) ,3 , d{d-W) ^2 , 31 + 22d-20d2+3d3 ^ , 2d-3^ , 
=5^aCa^Si^ +^^^ + 4(d-3)(d-5) ^ + + 

(3d - l)(d - 4) ^2 , 22 - 15d + 3^^ ^ ^ 3 + 



1 - 2d „ /d - 3 2 1 D M ^2 1 



+ + J'--{J+B.^^-£^-j^{JJ, (421) 

. j= - (1 «^ - 1 e . i) + i..^ e= + i..! , (4.22) 

+ i±i^4._.5+*5±i))B,, (4.24) 

H _ n^f d{2 + 33d - 12d^ + d^) ^ d(7 + 27d - ed^) ^3 



192 " 32 



d(ll + 29d-4d2) 2 2 - 21d - 33d2 + 8d3 Sd^X 

16 ^ + 4(d^5^ j^^' (^-^^^ 

„M_ 4, ^2/ ^(^-5) ,3 d(rf-9) ,2 I 2 + 13d-3d^ , d\ 

The three-loop diagrams depicted in Fig. (4) can be expressed as (diagram u ^ d^/Nc = 
dA(CA — 2Cf) contributes to two color structures) 

rf"^"' = ff^dAA^fCp^ (2(d - 1)(F2 - 2F3 + F5 + Fio) + 4F6 + (d - 5)^7) , (4.27) 
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Figure 4. Three-loop diagrams contributing to N[, and to Nf (last diagram only). 
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d- 5 
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Pif' = g''dKN,CA ( - 2Fi - - Fe + - Fg 



(4.28) 
(4.29) 
(4.30) 

(4.31) 



Summing all diagrams, the gauge parameter ^ explicitly drops out (and so do the 3-loop 
sum-integrals -B4 and S5), leaving the (bare) result 



The last three-loop diagram depicted in Fig. (4) gives 
p^} =g^dANf({d-f>)F^ + 



1 - d)dA + 4(2-"' - l)iVfC7A) lf+ 



+ g^dj—^ Ud - 1)Ca + 2(3 + 42-^^ - 24-'^)iVf) 



(1 - d)^ 



2((i - 5)^1 + ^2 + 2^3 + 



+ g^dANiCA 



l-d 



S{d - 5)Fi + 2{d - 3)Fe - {d - 5)Fr + 8^3 + 



,1 



l-d 



2{d - 1)(F2 - 2F3 + F5 + Fio) + 4F6 + (d - 5)F7 + 



+ (^^AiV/T 4((i - 5)F4 + {d- 3)F7 + 4F9 + 0(<7^) . 



(4.32) 



Note that also the spurious poles (at d = 3 — 2e) which are present in Eqs. (4.18), (4.19) 
and (4.21) (and which would have severely limited our ability to evaluate the constant term, 
requiring 3-loop sum-integrals to higher order in e) have canceled completely. Comparing 
with the corresponding expression given in Eq. (31) of [19], we note complete agreement, up 
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to two typos in that reference (in their 3-loop terms, Cp^-pX\l2L\ is missing a prefactor (1+e), 
and Cf7fXiTiX2 should be multiphed by 1/2; these correspond to our structures NiCp^F\ and 
A'^fCp-Fs above). 

Renormahzing the couphng in the MS scheme via 

2 -2e 2 2 /i /^O 5r , 4 \\ o HCa - 2iVf i a oo\ 

9=^ fbare = 9R(^l- — Y^ + ^(5r)J ' = § (^-33) 

from which (using that (7^are does not depend on the renormahzation scale /i) its running 

2 

fi'R(/^) ~ 5^,(^0) (l + 2/3o]|^ TT + ^(S'r)) follows, and expanding the master sum-integrals 
around d = 3 — 2e (see Appendix), the result coincides with the expression given in Eq. (32) 
of [19] (note that they write the MS scale parameter as A, where A^ = /x^ = 47re~'^^ /u^), 
which has subsequently been used to build pqcd(2^) according to Eq. (2.3). We will not 
re-iterate discussions about convergence as well as (renormahzation) scale dependence here, 
but instead only refer to the literature [18, 32-34]. Our main point was the re-derivation of 
Phard(2^) Eq. (4.32), demonstrating gauge-parameter independence and using automated 
computer- algebra methods, which are naturally extensible to higher orders. 

5 Conclusions 

We have re-derived the results of the original 3-loop computation [19] for the QCD pressure. 
While the original calculation was performed in Feynman gauge, we have used covariant gauge 
and shown explicit cancellation of the gauge parameter. As a slight improvement over [19], 
our IBP reduction revealed that the basis of non-trivial 3-loop master integrals used in the 
original computation could be reduced by one. 

In future work, we hope to be able to employ the IBP setup to the 4-loop level, thus con- 
tributing to the last missing piece needed for the physical leading- order (i.e. g^) determination 
of the pressure within the dimensionally reduced effective theory framework [24, 26, 32] - the 
only currently known framework allowing for a weak-coupling expansion of thermodynamic 
observables that is systematically improvable. 
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A One- and two-loop vacuum sum-integrals 

The one-loop bosonic tadpole is known analytically and reads (d = 3 — 2e) 

Tp (i^^)™ (27rT)2™-" 
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whereas the fermionic tadpole can be related to the corresponding bosonic one via 



-\2m—n—d 



(A.2) 



>{P} (P'r ' 

which immediately follows from the scaling relations explained in Sec. 3. 

As mentioned above, via integration-by-parts relations all two-loop integrals are express- 
ible in terms of products of two one-loop tadpoles which means they are also available ana- 
lytically up to arbitrary order in e. A typical reduction, following from IBP and being used 
in the 2-loop contribution to Phard 

(T), reads 

1 



S = 



This remarkable result in fact follows from the IBP relation 

1 



. 



(A.3) 







with fi = {d- 2){pi + q,) + — (Po + Qo)iqiPo - PiQo) , 



(A.4) 
(A.5) 



tpQ P^-"P^QHP-Qr 

2 

from which, after working out the derivatives, using the shift Q ^ P — Q and exploiting the 
P ^ Q symmetry, it follows that = {d — 2){d — 3)5, proving Eq. (A.3). Note that in the 
literature, by explicit integration Eq. (A.3) was only known to hold to 0{e) [18]. 



B Three-loop vacuum sum-integrals 

All non-trivial three-loop master sum-integrals that are needed for phard of Eq. (4.32) have 
been evaluated in Ref. [18], and been subsequently summarized in the literature [19, 22]. In 
the notation of Sec. 4, they read 

rpi / ,,2 \ 3e 



B2 

Fs 
Fg 
Fio 



1 1 



167r2 V47rT2 J e 24 



1 + (i - 37e + 8Zi - 2Z3) e + 0{e) 



167r2 V47rr2 ) 7 216 



l+(i-fT7E + ff^i-i?^3)6 + 0(e)^ 



rp4 



2 \3e 



1 1 



167r2 \A7rT^ J e 96 



1 + - 37e - f ln2 + 8^1 - 2Z3) e + 0{e)' 



2 \ 3e 



1 -29 



167r2 V47rr2 J e 1728 L 



l + (i-i7E-iln2 + ^Zi-i§Z3)6 + 0(e) 



167r2 V47rr2 ) 7 108 



35 _ 3 _ 63 
8 2 'B 10 



ln2 + 5Zi - iZ3)e + 0(e) 



2 \3e 



1 -1 r 
167r2 V47rT2 J 7 192 



1 + (If + 37e + f ln2 - 4Zi + 4Z3) e + 0(e) 



(B.l) 
(B.2) 
(B.3) 
(B.4) 
(B.5) 
(B.6) 



where Zj = yjzz^- Fq can be obtained via Eq. (3.8), and coincides with the expansion given 
in [19]. 
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